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Abstract 

Objective: To investigate the optimal control strategy for Plasmodium vivax 
malaria transmission in Korea. 

Methods: A Plasmodium vivax malaria transmission model with optimal control 
terms using a deterministic system of differential equations is presented, and 
analyzed mathematically and numerically. 

Results: If the cost of reducing the reproduction rate of the mosquito population 
is more than that of prevention measures to minimize mosquito-human contacts, 
the control of mosquito-human contacts needs to be taken for a longer time, 
comparing the other situations. More knowledge about the actual effectiveness 
and costs of control intervention measures would give more realistic control 
strategies. 

Conclusion: Mathematical model and numerical simulations suggest that the use 
of mosquito-reduction strategies is more effective than personal protection in 
some cases but not always. 



1. Introduction 

Malaria is a mosquito-borne infectious disease 
caused by a eukaryotic protist of the genus Plasmodium. 
Malaria is naturally transmitted by the bite of a female 
Anopheles mosquito. The primary vector in Korea is 
reported to be A sinensis. Since the re-emergence of 
Plasmodium vivax malaria in 1993 [1,2], it has been 



endemic and continues to cause extensive morbidity in 
Korea, despite the huge efforts invested to control it. 

The first mathematical malaria model proposed by 
Ross [3], was subsequently modified by MacDonald, 
which has influenced both the modeling and the appli- 
cation of control strategies to malaria [4]. Recently, the 
optimal control theory has been applied to malaria 
Okosun et al [5], and to vector-borne disease Lashari 



"Corresponding author. 
E-mail: yongkuk@knu.ac.kr 

This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License 
(http://creativecommons.Org/licenses/by-nc/3.0) which permits unrestricted non-commercial use, distribution, and reproduction in 
any medium, provided the original work is properly cited. 



Copyright © 2012 Korea Centers for Disease Control and Prevention. Published by Elsevier Korea LLC. All rights reserved. 



Optimal control strategy of Plasmodium vivax 



129 



Tnhlp 1 


\^r\f* Hf i cr > t"ii'\tii^Ti r\T nQrompfprc Tinv trip* mnnpl 
1 11 C UCaUl 1LH1U11 Ul Udl alllC IC1 S 1U1 111C 111UUC1 


Parameter 


Description 


°m 


J^CI ^dUlLa idlC Ul 11CW1V LllJLimiJL; tlUUll lllUaU UlLUCf> 


8 , 


Infected mosquito to human transmission efficiency 


Phm 


Infected human to mosquito transmission efficiency 


ff 


Average number of contact made to host by a single mosquito 


r 


Per capita rate of progression of humans from the infectious state to the susceptible state 


P 


Probability of exposed humans going through short-term incubation periods 




Per capita rate of progression of humans from the short term of exposed state to the infectious state 


± h 


Per capita rate of progression of humans from the long term of exposed state to the infectious state 



et al [6], who modified the model of Blayneh et al [7], 
but introduced some awkward terms. 

Models for Plasmodium flacipamm malaria or 
vector-borne diseases have been studied by many 
researchers [8—10]. In contrast, models for P vivax 
malaria are rare. Recently, Nah et al [11] proposed 
a model of P vivax malaria transmission. In this paper, 
by combining the ideas of Blayneh et al [7] and Nah 
et al [1 1], we propose the deterministic model of P vivax 
malaria transmission with optimal control terms. Using 
the optimal control theory, we sought optimal control 
strategies of P vivax malaria transmission in Korea. 

2. Materials and Methods 

2.1. Model description: optimal control 

To construct a deterministic model for P vivax 
malaria transmission with control terms, the model of 
Nah et al [11] was modified and optimal control terms 
inspired by the model of Blayneh et al [7] were added as 
follows: 



long term exposed (E H ), and infectious (///). Mosquito 
population M(t) is also divided into two classes: 
susceptible (Sm), and infectious (Im)- Note that the 
mosquito population M{t) is not constant while human 
population H(t) is constant. 

The factor of 1 — u\(t) reduces the reproduction rate 
of the mosquito population. It is assumed that the 
mortality rate of mosquitoes (susceptible and infected) 
increases at a rate proportional to u\(t), where p > 0 is 
a rate constant. In the human population, the associated 
force of infection is reduced by a factor of 1 — U2{t), 
where U2(t) measures the level of successful prevention 
efforts. In fact, the control U2(t) represents the use of 
prevention measures to minimize mosquito-human 
contacts. Table 1 lists detailed descriptions of the 
parameters. The system (1) has a unique solution set. 
(See Appendix A for detail.) 

An optimal control problem can now be formulated 
for the transmission dynamics of P vivax malaria 
transmission in Korea. The goal is to show that it is 
possible to implement time dependent anti-malaria 



' ^= - $ mh a{\ - u 2 (t))^I M (t) + rI H (t) 
-!L= p $ mh o(\ - U2 (t))^>I M (t) - r h E° H (t) 
-^=(1 -p)M\ - u 2 (t))^I M (t) - T{E' H {t) 
d -^=r h E s H {t) + T' h E> H {t)- r l H {t) 

-Jf =b r»(^ - ui{t)){S M {t) +/jf(0) _ PwC 1 - u 2(t))-^-Su(t) -b m S M {t) - pllx{t)S M {t) 
^ ^=/5 ta ff(l - Ui{t)P^-S M {f) - bJ M {t) - P u x {t)I M {t) 



In the model, human population H(t) is divided into control techniques while minimizing the cost of imple- 
four classes: susceptible (Sh), short term exposed (E S H ), mentation of such control measures. 
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Table 2. The parameter values for the model 



Parameter Value 

b m 0.7949 [0.1,1.5] 



$mh 


0.5 




0.5 


a 


0.3 [0.25,0.5] 


r 


0.07 [0.01,0.5] 


P 


0.25 


n 


1/25.9 


n 


1/360.3 



An optimal control problem with the objective cost 
functional can be given by 

T 

/B B 
{AI H {t)+^-u](t)+^-u\(t))dt, (2) 

o 

subject to the state system given by (1). 

The goal is to minimize the infected human pop- 
ulations and the cost of implementing the control. In the 
objective cost functional, the quantities A, B\ and B2 
represent the weight constants of infected human, for 
mosquito control and prevention of mosquito-human 
contacts, respectively. The costs associated with 
mosquito control and prevention of mosquito-human 
contacts are described in the terms B\u\ and B2U2, 
respectively. 

Optimal control functions (h*,m|) need to be found 
such that 

j(ll" x ,1*2) =XX\\X\{J{U\,U2)\(U\,U2)£ U}, 

subject to the system of equations given by (1), where 
U = {(ui,u 2 )\iii(t) ispiecewise continuous on [0, T],0 
<Ui(t) < l,i=l, 2} is the control set. 

Such optimal control functions (u* , uV) exist, and the 
optimality system can be derived. (See Appendix B for 
detail.) 



2.2. Numerical simulation 

Using the forward-backward sweep method, the 
optimality system was solved numerically. This consists 
of 12 ordinary differential equations from the state and 
adjoint equations, coupled with the two controls. In 
choosing upper bounds for the controls, it was supposed 
that the two controls would not be 100% effective, so 
the upper bounds of u\ and 112 were chosen to be 0.8. 
The weight in the objective functional is A\ = 1000. The 
parameters in Table 2 were adopted from other articles 
[11] and used for our simulation. 

We simulate the model in different scenarios. Figure 1 
depicts scenarios for the state variables of the model for 
the case when the cost is the same for the two controls. 
Figure 2 depicts scenarios for the state variables of the 
model for the case when the cost of prevention measures 
to minimize mosquito-human contact is more expensive 
than the cost of reducing the reproduction rate of the 
mosquito population. Figure 3 depicts scenarios for the 
state variables of the model for the case when the cost of 
reducing the reproduction rate of the mosquito population 
are more expensive than the cost of prevention measures 
to minimize mosquito-human contacts. 

It is also worth noting that different initial mosqui- 
toes populations do not have effect on the optimal 
strategies (Figures 4 — 6). 

2.3. Results 

If the cost of reducing the reproduction rate of the 
mosquito population is more than that of prevention 
measures to minimize mosquito-human contacts, the 112 
control needs to be taken for a longer time, comparing 
the other situations (Figures 1 to 3). In that situation, full 
effort for 112 is needed after the high peak of infected 
human population. 

On the other hand, Figures 4 to 6 suggest that even 
though the mosquito population is not so high in initial 
point, full efforts for u\ and 1/2 are needed for at least 
some of the time. 
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Figure 4. Optimal controls when B\=Bi = 1000 with low mosquito population. 
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Figure 5. Optimal controls when Bi = W, .82 = 1000 with low mosquito population. 
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Figure 6. Optimal controls when i?i = 1000, B-2 = 10 with low mosquito population. 



3. Discussion and Conclusions 

After 1993's reemergence of malaria, the endemicity 
of P vivax malaria is becoming a growing concern in 
South Korea. Public health advisories were subsequently 
issued to apply community mosquito control and 
personal protection. 

The purpose of this work is to suggest optimal control 
strategies of P vivax malaria in different scenarios. In all 
cases, optimal control programs lead effectively reduce 
the number of infectious individuals. We have used 
a deterministic model with time-dependent parameters to 
develop the transmission dynamics of P vivax malaria in 
Korea. For numerical simulations, most parameters were 
adopted from other articles [11]. 

Mathematical model and numerical simulations 
suggest that the use of mosquito-reduction strategies is 
more effective than personal protection in some cases 



but not always. Public health authorities should choose 
the proper control strategy where their situation lies in 
the scenarios discussed in the Results section. 
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Appendix A. The existence and uniqueness 
of solution 

We consider system (1). We obtain the existence and 
uniqueness of solution. In here we are given a suitable 
control set. 
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Theorem 1. The system (1) with any initial condition 

has a unique solution. 

Proof. We can rewrite (1) as : 



dX_ 
dt 



■■AX + F(X, U), 



where X= [S H ,E S H ,E' H ,I H ,S M ,I M ] T , 





0 


0 


r 


0 


0 \ 


0 


1 h 


0 


0 


0 


0 


0 


0 


'T'S 


0 


0 


0 


0 


n 


n 


— r 


0 


0 


0 


0 


0 


0 


-b m 


0 


Vo 


0 


0 


0 


0 


-bj 



U=[u u u 2 ] T and F(X,U)-- 



-Pmh <7 ( l-"2 (0 l M W > 

^ ff (l- M2 (0^ M*),(l-P)/W (l-"2«^ 

4? (0 ,0 /> m ( l-M , (?) )(% (0 Vhm°{ l-»2 (f) ) 



\\F(X U U)-F{X 2 ,V)\\ 

< fi mh a/H\{\ - u 2 )S H \Im\ - (1 - U2)5Wm2| 
+pP m ho/H\{\ - u 2 )ShJmi ~ (1 - i^S/h^I 

+(1 -p^a/HlH -u 2 )S m I Ml - (l-v 2 )S H2 I M2 \ 
+b m \(l -u 1 )(S M i +Imi) ~ (1 -i>i)(Sm2 + -Jm2)| 
+/? Am <r///'|(l - u 2 )S M Jm - (1 - v 2 )S m %Ih2 \ 
+p\u\Sm -viS M2 \+fihm^/H\{\ - u 2 )S m Jh\ 

-(1 -V 2 )S«2.Zh2|+P|M1-Jm1 -Vl/jtf 2 | 

= 2^/^1(1 - K 2 )S m / K i - (1 - v 2 )S H2 I M2 \ 
+b m \(l - ui)(S M \ +Imi) - (1 -Vi)(S M 2+Im2)\ 
+2(5 hm o/H\(\ - u 2 )S M \Im - (1 - v 2 )S M 24d 
+p|wi5 M i - viS M2 \ +p\uJmi - vJmz\ 
We estimate the 4 terms in the right side of (i): 

2/? mA cr///|(l - u 2 )S m I Ml - (1 - i^S/fiAd 

< 2ft rt ff/J?[2M'(0)|Sin - 5ml + 2//|/mi - Im\ 
+HM(0)\u 2 -v 2 \], 

b m \(l -ui)(Sm+Im) ~ 0 - ui)(Sm 2 + 
<2fe m [|5 M1 -S M2 \ + \I M1 -/jfll+AfCOJIiii-Wjl], 

2f} hm o/H\{\ -u 2 )S H \Im\ - (1 —v 2 )SmIm\ 

< 2fS hm a/H[2M(0) \I m -I H2 \ + 2H\S m 
+M{0)H\u 2 -v 2 \], 



(i) 



H 

1 T 



H 



4r(0 



. So let G(X,(7)=A¥+F(X,C/). Defined matrix pMj#i - </»E^(0)|«i -Vi| + 



^ is a linear. So A is a bounded operator. Define 
a matrix norm and a vector norm as follows 
I 1 1 = 2,- 7 - |«,y I ,| I A- 1 1 = 1 1 (x,-) 1 1 = 2,- 1 , respectively. To 
show the existence of solution of the system (1), we 
have to prove that F(X,U) satisfy a Lipschitz condition. 
Let 

H(t) : =S„(t) +E s H (t) +E' H (t) +I H (t). 
and 

M(t) : =S M (t)+I M (t). 
But 

j t H'(t) = 0. So, H(t)=H (a constant^ < 00 , 



(11) 



(iii) 



(iv) 



(v) 



and 



dt 



M(t) = - ui (/) (b m + p)M(t) < 0. Hence, 

M(t) < M(0) < 00. 
For any given pairs (X\ ,U),(X 2 ,V), 



s hj , E s Hj , E H . , I HJ , S M j , Imj 



,7=1,2, 



p\uJ M \ -vil m \ < p[M(0)\ui -vi\ + \I m -ImiW- (vi) 
Since 0 < m,,u, < 1,0 < ShuIhi < H and 0 < Sm^m < 
M(0),z = 1,2, by (i)-(vi), we have 

||F(*i,C/)-F(Jr 2> r)|| 

<(4/? m , ff M(0) ////) \S m -S H2 \+(4fi hm aM(Q) ///) 1^, -I H2 \ 
+{2b m +4fi hm aM{Q)+p) \S m -S M2 \ 
+(2b m +4p mh crM{Q)+p)\I Ml -I M2 \ 
+{2b m M(0)+2pM(0))\u l -v l \ 

+ (2p mh +2!3 hm aM(0))\u 2 -v 2 \ 

<K{\Sm Sh2 \-\-\Sm 1 — Sm2 — Im |+|/mi — ^M2 1 

+|m 1 -u 1 |+|m 2 -u 2 |) 

<^(||^-A 2 ||+||f/-F||), 

where K = maz{4P mh (rM(0)/H,4P hm oM(0)/H,2b m 
+4(3 hm trM(O) + p, 4/3 mA ff +2fe m + 
p,2M(0)(A m + p),2 l 8 m/! +2 l 8 /im o-M(O)} Thus, we obtain 
F(X,U) is uniformly Lipschitz continuous. Let 
Z=||yf|| < 00. Then, 



U = (u u u 2 ) T , V = (v u v 2 ) T , 
we obtain, 



HGCYi, t/) - G(X 2 , V)\\ < (K+EjQXi -X 2 \ 

+ \\u-v\\). 
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Hence, the system (1) satisfy all conditions of the 
Picard-Lindelof Theorem ([12,13]) and also the function 
F(X, U) is continuously differentiable. Therefore, the 
system (1) have a unique solution. 



Appendix B. Analysis of optimal control 
control problem 

We are to prove the existence of optimal control pairs 
for the system (1). Firstly, We set control space 

U = { (u\ , u 2 )\ui is piecewise continuous on [0, T] , 0 
<u t (t) < 1,1=1,2}. 

We consider an optimal control problem to minimize 
the objective functional: 

T 

J{u u u 2 ) = J (AI H {t) > Bl -*<* > Bl ~ 2 < 



t)+^u]{t)+—u\{t) )dt. 



Theorem 2. There exist an optimal control u\ and u\ 
such that 



j ( hi | , liry 



min J{u\,u 2 ) 



,M 2 E(7 



(3) 



subject to the control system (1) with initial conditions. - 
Proof. To prove the existence of an optimal control pairs 
we use the result in [14]. The set of control and corre- 
sponding state variables is a nonempty. Because for 
each control pairs we have proved in the Theorem 1 that 
there exists corresponding state solutions. And also it is 
ok when the control u\=ui = §. Note that the control 
and the state variables are nonnegative values. The 
control space U is close and convex by definition. In the 
minimization problem, the convexity of the objective 
functional in u\ and Uj have to satisfy. The integrand in 



B\ 7 

the functional, Aljj(t) + — u\ (t) 



B 2 2 



u 2 (t) is convex 



function on the control u\ and 112 ■ Also we can easily 
check that there exist a constant p > 1, a numbers cj\ > 
0 and a>2 > 0 such that 

J(u U U 2 ) > tol + W2(|wi| 2 + |"2| 2 ) P/ 

which completes the existence of an optimal control. To 
find the optimal solution we apply Pontryagin's 
Maximum Principle ([15—17]) to the constrained 
control problem, then the principle converts (1), (2) and 
(3) in to a problem of minimizing pointwise a Hamilto- 
nian, Jf , with respect to u\ and The Hamiltonian for 
our problem is the integrand of the objective functional 
coupled with the six right hand sides of the state 
equations: 



B 



JT(X(0, u(0, A(0) =AI„(t) + ^ u\ (t) + -fu z 2 (t) 



B 



2 . .2/ 



(4) 



where gi is the right hand side of the differential equa- 
tion of the rth state variable and also 



X(f) = (Sh,E s h ,E 1 h ,I h ,S m ,Im) , U(0 = (hi (t) , « 2 (0) an d 
A(f) = (Ai (?) , A 2 (0 , A 3 (t) , A 4 (?) , A 5 (0 , A 6 (0) • 
By applying Pontryagin's Maximum Principle([18]) if 
(x*(t),U*(t)) is an optimal control, then there exists 
a non-trivial vector function X(t) satisfying the 
following equalities: 

dx_dJf(x(t),u(t),A(t)) 
dt 



0 = 



9A 

3jr(x*(f),u*(0,A(0) 



du 

3^(x«(0,u*(0,A(Q) 
9x 



If follows from the derivation above 



.3^ 
M r = o,if— <o 

oui 

0<u*< l,if— =0 
.„9<Jf 

M*=l,lf— >0. 



Now, we apply the necessary conditions to the Hamil- 
tonian 3^. 



Theorem 3. Let S* H (t),E%(t),E l *(t),I* H (t),S* M (t) and 
I^i(t) be optimal state solutions with associated optimal 
control variables u\ and u% for the optimal control 
problem (1) and (2). Then, there exist adjoint variables 
h(t), fa{t), A 3 (f),A4(f),A 5 (f) and A 6 (?) that satisfy 



r i (t) = X 1 (t)B mh( r(l-u*(t)) II 



1 

H 



-X 2 {t)pP mh a{\~u* 2 {t)) II 



H 



H 



-h(t)(l-p)P mh a(l-u* 2 {t)) I M 

Ut)=h{t)n-x,{t)r h 



X 3 (t)=X 3 (t)T l h -h(t)T> h 



X' 4 (t)=-A-X 1 (t)r+X 4 (t)r+X 5 (t)8 hm a(l-t4(t))Sl 



H 



-h{t)fa m (r{\-u* 2 {i))S* M 



H 



X' 5 (t)= -X 5 (t)(b m (l-u*(t)) -f3 hm a(\-u* 2 (t)) I* H 
-b„,-pu*(t)) 



H 



-X 6 {t)R hm a(\-u* 2 {t)) F H — 
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X' 6 (t) = X 1 (t)P ml ,a(\-u* 2 (t))S* H - 



H 



-X 2 {t)p$ mh a{\-u* 2 {t)) S* H 



H 



-X 3 (t)(l-p)P mh cr(l-u* 2 (t))S* H - 

-X s (t)b m {\ -u\(t)) + X 6 (b m + puf) 
with transversality conditions(or boundary conditions) 

Xj(T) =0,7 = 1,2,..., 6. (5) 
Furthermore, the optimal control u\ and u% are given by 



u\{t) =min|max|o,-^- ((p + b m )X 5 (t)S* M + (pX 6 (t) + b m X 5 (t))I^(t)) J, 1 

M *(*)=mmjmaxjo,-^ ^( - A,(?) + P X 2 (t) + (1 -p)X 3 (t))/3 ml ,aS , H (t)r M (t)^ 

1 



(6) 



Proof. To determine the adjoint equations and the Using the property of the control space, we obtain the 

transversality conditions, we use the Hamiltonian (4). characterizations of u\{t) and m^M m F rom 

From setting SH{t)=S* H {t),E s H {t)=E s ^{t),E l H {t)= the fixed of start time, we have transvesality 

E'i{t), I H (t)=I* H (t),S M (t)=S* M (t) and I M (t)=I* M (t), conditions (5). 
and also differentiating the Hamiltonian (4) with 
respect to Sh,E s h ,E 1 h ,Ih,Sm and /M,we obtain 



a'. 



X' 



X' 



dffl 11 1 

— =X l /3 mh a(l -u* 2 )r M --X 2 pP mh a(\ -u* 2 )l* M — - A 3 (l - p)P mh a(\ - u* 2 )l* M — 

9,^ 



dE' H 

9,?r 



■■x 2 T s h - x 4 r h 



Xt,T 1 , — XaTI 



1 3 h 4 fi 



K = ~-^j- = ~ A ~ hr + X*r + X 5 $ hm o(\ -u* 2 )S* M -^- X 6 (3 hm a(\ -u* 2 )S M -^ 



XL 



= A 5 fb m (l -«*) ~fS hm a(l -u*)l*~ 



pu x 



1 

H 



X' 



9,#' 
ds M 

dffl I I 1 

XiP mh (r(l -u* 2 )S* H — - X 2 pP mh a(\ -u* 2 )S* H — - A 3 (l - p)P mll a(l - u* 2 )S H — 



9/, 



-A 5 6„,(l -«*) + A 6 (6„, + 
By the optimality conditions, we have 



o=^=B l u\(t)-x 5 (t)(b l „(s t M +i^ + ps^)-x 6 (t)pr M (t) 

d ffl 1 1 

0 = ^ U2 =B ^ +^(0^<r5*(/)/;(/)--A 2 (l)A, k <r5*W/;(/)- 



136 



B.N. Kim, et al 



References 

1. Feighner BH, Pak SI, Novakoski WL, et al. Reemergence of 
Plasmodium vivax malaria in the Republic of Korea. Emerg Infect 
Dis 1998 Apr-Jun;4(2):295-7. 

2. Ree HI. Studies on Anopheles sinensis, the vector species of vivax 
malaria in Korea. Korean J Parasitol 2005 Sep;43(3):75— 92. 

3. Ross R. The prevention of malaria. 2nd ed. London: Murray; 1911. 

4. MacDonald G. The epidemiology and control of malaria. London: 
Oxford University Press; 1957. 

5 . Okosun KO, Ouifki R Marcus N. Optimal control analysis of a malaria 
disease transmission model that includes treatment and vaccination 
with waning immunity. Biosystems 2011 Nov;106(2-3):136— 45. 

6. Lashari AA, Zaman G. Optimal control of a vector borne disease 
with horizontal transmission. Nonlinear Anal: Real World Appl 
2012 Feb;13(l):203-12. 

7. Blayneh KW, Gumel AB, Lenhart S, et al. Backward bifurcation 
and optimal control in transmission dynamics of West Nile virus. 
Bull Math Biol 2010 May; 72(4): 1006-28. 

8. Blayneh KW, Cao Y, Kwon H-D. Optimal control of vector-borne 
diseases: treatment and prevention. Discrete Contin Dyn Syst B 
2009;11:587-611. 

9. Bowman C, Gumel AB, van den Driessche P, et al. A mathe- 
matical model for assessing control strategies against West Nile 
virus. Bull Math Biol 2005 Sep;67(5):l 107-33. 



10. Pongsumpun P, Tang IM. Mathematical model for the trans- 
mission of P. falciparum and P. vivax malaria along the 
Thai-Myanmar border. Int J Biol Life Sci 2007;3(3):200-7. 
Summer. 

1 1 . Nah K, Kim Y, Lee JM. The dilution effect of the domestic animal 
population on the transmission of P. vivax malaria. J Theor Biol 
2010 Sep;266(2):299-306. 

12. Hartman P. Ordinary differential equations. New York: John 
Wiley & Sons; 1964. 

13. Kelley E, Petterson A. The theory of differential equations: clas- 
sical and qualitative. Pearson Education, Inc.; 2004. 

14. Lukes DL. Differential equations: classical to controlled. In: 
Mathematics in science engineering, vol. 162, New York: 
Academic Press; 1982. 

15. Kamien MI, Schwartz NL. Dynamic optimization: the calculus of 
variations and optimal control in economics and management. 
Amsterdam: North-Holland; 1991. 

16. Lenhart S, Workman JT. Optimal control applied to biological 
models. In: Mathematical and computational biology series. 
London/Boca Raton: Chapman & Hall/CRC Press; 2007. 

17. Sethi SP, Thompson GL. Optimal control theory - applications to 
management science and economics. Springer; 2000. 

18. Pontryagin LS, Boltyanskii VG, Gamkrelidge RV, et al. The 
mathematical theory of optimal processes. New York: John Wiley 
and Sons; 1962. 



